Extreme Value Statistics of the Total Energy in an 
Intermediate Complexity Model of the 

Mid-latitude Atmospheric Jet. 
Part II: trend detection and assessment. 

Mara Felici^'^, Valerio Lucarini\ Antonio Speranza\ Renato Vitolo^' * 

February 2, 2008 

^ PASEF - Physics and Applied Statistics of Earth Fluids, 
Dipartimento di Matematica ed Informatica, Universitd di Camerino 

^ Dipartimento di Matematica U. Dini, Universitd di Firenze 



* Corresponding author address: Dr. Renato Vitolo, Department of Mathematics and Informatics, Uni- 
versity of Camerino, via Madonna delle Carceri, 62032 Camerino (MC), Italy. 
E-mail: renato.vitolo@unicam.it 



1 



Abstract 

A baroclinic model for the atmospheric jet at middle-latitudes is used 
as a stochastic generator of non-stationary time series of the total energy of 
the system. A linear time trend is imposed on the parameter Te, descriptive 
of the forced equator-to-pole temperature gradient and responsible for set- 
ting the average baroclinicity in the model. The focus lies on establishing 
a theoretically sound framework for the detection and assessment of trend 
at extreme values of the generated time series. This problem is dealt with 
by fitting time-dependent Generalized Extreme Value (GEV) models to se- 
quences of yearly maxima of the total energy. A family of GEV models is 
used in which the location jj, and scale parameters a depend quadratically 
and linearly on time, respectively, while the shape parameter is kept con- 
stant. From this family, a model is selected by using diagnostic graphical 
tools, such as probability and quantile plots, and by means of the likelihood 
ratio test. The inferred location and scale parameters are found to depend in 
a rather smooth way on time and, therefore, on Te- In particular, power-law 
dependences of fi and a on Te are obtained, in analogy with the results of 
a previous work where the same baroclinic model was run with fixed values 
of Te spanning the same range as in this case. It is emphasized under which 
conditions the adopted approach is valid. 

PACS: 02.50.Tt, 02.70.-c, 47.1 l.-j, 92.60.Bh, 92.70.Gt 
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1. Introduction 



In the context of Climate Change, an intensely debated question is whether the statistics of 
extreme meteo-climatic events is changing (and/or will change) and, in case, how fast it is 
changing (and/or will change). For example, the role of time-dependence in the statistics 
of extreme weather events has been at the heart of discussions about climate change since 



the work by lKatz and Brownl (119921) . In particular, the detection of trends in the frequency 



of intense pr ecipitation has be e n the object of much r esearch, particula rly at regional 



(12002 . 



evel, see e.g . iKarl et all (119961) : Karl and Knight! (Il998h for the USA and lBrunetti et al. 



2004|) for the Mediterranean area. The general relevance of the problem has been 



highlighted in the 2002 release of a specific IPCC report on Changes in extreme weather 
and climate events (available at |http ://www.ipcc . ch/pub/support .htm) . In fact, the empha- 



sis laid on the subject by the IPCC report reverberated in many countries the question "is 
the probability of major impact weather increasing?". This question reached the big pub- 
lic almost everywhere and innumerable studies of trends in series of "extremes" were un- 
dertaken. These studies mainly deal with variables of local character, typically precipita- 
tion and temperature at specific stations. Moreover, most studies a re regional: see e.g. the 



proceedings of the Italia-USA meeting held in Bologna in 2004 (|Diaz and Nannill2006h 



for the relevance of the extreme events in the Mediterranean Climates and the INTERREG 
IIIB - CADSES project HYDROCARE, http://www.hydrocare-cadses.net, 
for impacts of extreme events in the hydrological cycle of the central-eastern Europe. 
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In a preceding, companion paper (IFelici et alj2006l) (which we refer to by Part I in the 



sequel) we have addressed the problem of extreme value statistical inference on statisti- 
cally stationary time series produced by a dynamical system providing a minimal model 
for the dynamics of the mid-latitudes baroclinic jet. There reported is, from mathematical 
literature, a suitable, rigorous, "universal" setting for the analy sis of the extrem e events 



in stationary time series. This is based on Gnedenko's theorem (|Gnedenkdl 1 9431) accord 



ing to which the distribution of the block-maxima of a sample of independent identically 
distributed variables converges, under fairly mild assumptions, to a member of a three- 
real paramet er farnily of distributions, the so-called Generalized Extreme Value (GEV) 



distribution (|Colesll200lb . The GEV approach to the analysis of extremse requires that 



three basic conditions are met, namely the independence of the selected extreme values, 
the consideration of a sufficiently large number of extremes, the selection of values that 
are genuinely extreme. This could be performed relatively easily for the case at hand. 

Part I was originally motivated by the interest in weather having "major impact" (on 
human life and property) in the Mediter ranean area, in particular in t ense precipitation and 



heat w a ves over Italy. See, for example. iBrunetti et al.l(l2002ll2004l):lLucarini et al 



2006bh : ISperanza et al.l(l2006h : ISperanza and Tartaglionel(l2006h : lTartaglione et all (120061) 



(2004, 



and the MEDEX Phase 1 report (available at http : / /medex . inm . uib . es /) for re- 
lated results and activities. The study reported in Part I has revealed, among other things, 
that diagnostics of extreme statistics can highlight interesting dynamical properties of 



the analyzed system. Properties which, thanks to the "universality of the GEV", can be 



investigated in a low dimensionality space of parametric probability density functions, 
although at the expenses of the total length of the observational record of the system in 
order to capture a sufficient number of independent extremes. A key role (that is presently 
being explored elsewhere, in the context of general atmospheric circulation theory) was 
played in Part I by the smoothness of variation of the extreme statistics parameters (aver- 
age, variance, shape factor) upon the external (forcing) parameters of the system. In this 
paper, again, we devote attention to exploring the statistics of extremes as a dynamical 
indicator, this time in the framework of the (typically meteorological) statistical inference 
problem of detecting trends in observations. 

The definition of a rigorous approach to the study of extremes is much harder when 
the property of stationarity does not hold. One basic reason is that there exists no univer- 
sal theory of extreme values (such as e.g. a generalization of Gnedenko's theorem) for 
non-stationary stochastic processes. Moreover, in the analysis of observed or syntheti- 
cally generated sequences of data of finite length, practical issues, such as the possibility 
of unambiguously choosing the time scales which defines the statistical properties and 
their changes, become of critical importance. Nevertheless, GEV-based statistical mod- 
eling offers a practical unified framework also for the study of extremes in non- stationary 
time series. In the applications, the three parameters of the GEV distribution are taken 
as time -dependent a nd time is introduced as a covariate in the statistical inference pro- 



cedure (Coles 



20011). The practical meaning of this assumption is that the probability of 



occurrence (chance) of the considered extreme events evolves in time pretty much as we 



are inclined to think in our everyday life. However, giving a scientific meaning to such 
an assumption is possible only in an intuitive, heuristical fashion: in an "adiabatic" limit 
of infinitely slow trends (but rigorously not even in such a limit). We adopt this point of 
view not only because it is in line with the common practice and view of extremes, but 
also because interesting dynamical properties can be inferred from extremes, in analogy 
with the findings in Part I. 

In the present paper we perform and assess time-dependent GEV inference on non- 
stationary time series E{t) of the total energy obtained by the same simplified quasi- 
geostrophic model that was used in Part I. The model undergoes baroclinic forcing to- 
wards a given latitudinal temp erature profile controlled by the forced equator-to-pole 



temperature difference Te', see 



Lucarini et al. 



( 2006 J 



dl) for a thorough description. We 



analyze how the parameters of the GEV change with time when a linear trend is im- 
posed on the large scale macroscopic forcing Te, that is, when Te is taken as a (linear) 
function of time. Since this functional relation is invertible, we derive a parametrization 
relating the changes in the GEV to variations in Te (instead of time). One major goal 
here is to present a methodological framework to be adopted with more complex models 
and with data coming from observations, as well as an assessment of the performance 
of the GEV approach for the analysis of trends in extremes in the somewhat grey area 
of non- stationary tirn e series. Methodologically, our set-up is somewhat similar to that 



of 



Zhang et al.l (|2003j) regarding the procedures of statistical inference. However, in this 



case we face two additional problems: 



1. as in Part I, the statistical properties of the time series E(t) cannot be selected 
a priori: in the stationary case {Te constant in time) and much less in the non- 
stationary case there is no explicit formula for the probability distribution of the 
observable E{t); 

2. moreover, in the non-stationary case we even lack a definition (and in fact a mere 
candidate) of what might be the probability distribution of E(t): ce rtainly not a 



freque ncy limit for t oo (and not by construction, as opposed to 



Zhang et al. 



20031) who use genuinely stochastic generators). 



This also means that we have no hypothesis concerning the functional form of the trend in 
the statistics of extremes of E(t), resulting from the trend imposed on the control param- 
eter Te- The lack of a general GEV theorem for non-stationary sequences implies that 
the choice of the time-dependent GEV as a statistical model is, in principle, arbitrary: 
other models might be equally (or better) suitable. Here comes into play the "adiabatic- 
ity" hypothesis mentioned above, which leads to the central statement of this paper: if 
the trend is sufficiently slow and if the statistical behaviour of the atmospheric model has 
a sufficiently regular response with respect to variations in the external parameters, the 
GEV remains a suitable model for inference of trend in extremes. 

The structure of the paper follows. In Sec. [2] we describe the general problem of the 
characterization of statistical trends in deterministic models, with both its conceptual and 
practical implications. Then in Sec. [3] we describe how the GEV modeling can be applied 
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to non-stationary time-series and how the quality-check of the fits is peri'ormed. In section 
Sec. m we present the time series considered in this work and the set-up of the numerical 
experiments performed with the atmospheric model. The inferences for various values of 
the trend in the forcing parameter Te are presented in Sec. [5] and a sensitivity analysis is 
carried out in Sec. [6l Comparison with the stationary case analyzed in Part I is given in 
Sec. U\ In Sec. [H finally, we summarize the main findings of this work, highlighting the 
future research developments. 



2. Statistical trends: the theoretical problem 



The stochastic generator used in this paper to produce the time series is a deterministic 
model (an ordinary differential equation), whose dynamics, for the considered range of 
values of Tr, is chaotic in the sen se th at it takes place on a s trange attractor A in phase 



space (lEckmann and Ruelle 



1985b . See 



Lucarini et al 



(2006c 



d|) for a study of the proper- 



ties of this attractor, including sensitivity with respect to initial conditions. The statistical 
behaviour of this type of time series is determine d by the Sinai-Ruelle-Bowen (SRB) 



probability measure fj, (|Eckmann and Ruellelll985b : this is a Borel probability measure 
in phase space which is invariant under the flow /* of the differential equation, is er- 
godic, is singular with respect to the Lebesgue measure in phase s pace and its co nditional 



measures along unstable manifolds are absolutely continuous, see 



YoungI (120021) and ref- 



erences therein. Moreover, the SRB measure is also a physical measure: there is a set V 
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having full Lebesgue measure in a neighbourhood f/ of A such that for every continuous 
observable (j) -.U ^ M, we have, for every x eV , the frequency-limit characterization 



1 







lim- / 0(f(x))rft= / <pdii. (1) 

t — ^OO I 



Existence and uniqueness of an SRB measure ji have been proved only for very 
classe s of flows /* (in particular, for flows that possess an Axiom-A attractor, see 



special 



Young 



(|2002r) '). However, existence and uniqueness of /i are necessary to define a stationary 
stochastic process associated to an observable (j). In turn, this allows to consider a given 
time series of the form cjy^f^^x)) : t > as a realization of the stationary process, justi- 
fying statistical inference on a solid theoretical basis. In part I, we conjectured existence 
an uniqueness of an SRB measure for the atmospheric model, providing the theoretical 
foundation to the application of GEV models in the inference of extreme values. 

In certain cases of non-autonomous ordinary differential equations (for example, if the 
dependence on time is periodic), it still is possible to define, at least conceptually, what 
an SRB measure is. However, in the present case, due to the form of time-dependence 
adopted for the parameter Tg, the atmospheric model admits no invariant measure. This 
means that there is no (known) way to associate a stochastic process to the time series 
of the total energy. In other words, it is even in doubt what we mean by "statistical 
properties" of the time series, since it is impossible to define a probability distribution. 
This conceptual problem has a very serious practical consequence: the "operational" 
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definition of probability as a frequency limit (as in ([T])) is not valid in the non-stationary 
case, since the time evolution is not a sampling of a unique probability distribution. Even 
if one assumes the existence of a sequence of distinct probability distributions, one for 
each observation, one realization (the time series) does not contain sufficient statistical 
information, since each distribution is very undersampled (with only one observation). 

Despite all these problems, the results in Part I suggest a framework which is, for 
the moment, formulated in a heuristic way. Suppose you evolve an initial condition x in 
phase space by the flow /* of the autonomous atmospheric model, that is the system in 
which Te is kept fixed to some value T^. After an initial time span (transient), say for t 
larger than some to > 0, the evolution f*{x) may be thougth to take place on the attractor 
A and time-averages of the form 



which is the average of (p by the SRB measure /io existing at the value Te = (the 
"attractor average" at T^). Now suppose that at some ti > the parameter Te is abruptly 
changed to some value > T^: there will be some transient interval, call it [^1,^2], 




(2) 



may be considered as approximations of 




(3) 
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during which the evolution f^{x) approaches the new attractor of the atmospheric model, 
that is the attractor existing for Te fixed at T^. After that time span, the evolution may 
be considered to take place on the new attractor. 

In our case, though Te varies continuously (linearly) with time, if the trend magnitude 
is low, then Te may be considered constant (with good approximation) during time spans 
that are sufficiently long in order to have both convergence to the "new" attractor and good 
sampling of the "new" SRB measure, in the sense sketched above. Though admittedly 
heuristic, this scenario allows to clarify under which condition it is still reasonable to 
speak of "statistical properties" of a time series generated by a non-autonomous system: 
namely, the closeness to a stationary situation for time spans that are sufficiently long. 
This is the "adiabatic" hypothesis which we mentioned in the introduction. An essential 
ingredient for this to hold is that the statistical properties of the autonomous model do 
not sensibly depend on the external parameter Te, in the sense that no abrupt transitions 
(bifurcations) should take place as Te is varied. This was indeed checked for the system 
at hand in Part I. Notice that the validity of the "adiabatic" hypothesis also has a useful 
practical consequence: one can use the statistics of the stationary system as a reference 
against which the results for the non-stationary case can be assessed. Having this scenario 
in mind, we proceed to the description of the time-dependent GEV approach in the next 
section. 
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3. GEV modelling for non-stationary time series 

3a. Stationary case 



The GEV approach for sequences of independent, identica 



variables is by now rather standard dCasti. 



Falk et al 



1994 



Galambos 



Reiss and Thoma: 



1 



2001 



1978 



Gumbel 



Tiago de Qliveira 



lo 



1988 



1958 



rem (|Gnedenkdll943l : lFisher and Tippettill928h . Consider a time series, assumed to be a 



Coles 



ly di s tributed (i.i.d.) random 



2001 



Jenkinson 



1955 



Embrechts et al. 


1997 


; Lindgren et al. 


1983 



19841). The foundation is Gnedenko's theo- 



realization of a sequence of i.i.d. random variables. The time series is divided into m con- 
secutive time-frames (data blocks), each containing n subsequent observations, equally 
spaced in time. Denote hy zi, . . . , Zm the sequence of the block maxima taken over each 
time-frame. Under fairly mild assumptions, the distribution of the block-maxima con- 
verges, in a suitable limit involving a rescaling, to the GEV distribution, defined as 



G{x; fi, a, = exp 



1 + e 



X — jJ 

a 



(4) 



for all X in the set {x : 1 + ^{x — fj) / a > 0} and G{x) = otherwise, where a > and 
^ 7^ 0. If ^ = the limit distribution is the Gumbel distribution 



G{x; /i, cr, 0) = exp ( — exp 



X — jl 



a 



X G 



(5) 
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GEV inference consists in estimating the distributional parameters (/i, a, ^) (called lo- 
cation, scale and shape parameter, respectively) from the available data. A widely used 
technique consists in numerically maximizing the log-likelihood function a, For 
^ 7^ 0, this is defined as 



- m log cr - ( 1 + ^ <j log 



t=i 



1+e 



Zt- fJ' 

a 



1 + e 



a 



(6) 



while an analogous formula holds for ^ = (|Colesll200lh . 



A generalization of the GEV theorem holds for time series that are realizations of 
station ary stochastic processes such that the long-range dependence is weak at extreme 



levels dLeadbetter 



1974 



19831) . In the applications, this property is assumed to hold 
whenever the block maxima are uncorrected for sufficiently large block sizes. In this 
case, GEV inference and assessment is carried out by t he sam e tools (maximum likeli- 



hood, diagnostic plots, etc.) used in the i.i.d. context, see 



Coles 



( 2001 ). 



3b. Time -dependent case 

If stationarity of the time series does not hold, then the limiting distribution function is no 
longer bound to be the GEV or any other family (Hj): no theories of extreme values exist 
in this context. Some exact results are known only in certain very specialized types of 
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non-stationarity (|Husleij|1986tlLindgren et al.lll983|) . but it is very unlikely that a general 



theory can be established. However, GEV-based statistical modeling of extreme values 
can be performed also in the case of time-dependent phenomena by adopting a prag- 
matic approach, where the GEV distribution (HI) is used as a template: time-dependent 
parameters and a(t) are considered, yielding a GEV model of the form 



G(x;/i(t),cr(t),0. 



(7) 



Usually ^ is kept time-independent in order to avoid nume rical problems, since it is 



20011) . Different kinds of time- 



usually the most delicate parameter to estimate (iColesI 
dependence can be imposed for cr{t) ). In this paper, we adopt a simple polynomial 
family of models: 



= fJ,0 + /ilt + /i2t^, Cr(t) = (To + ait, 



(8) 



with /io,i,2 and o-q,! G R. GEV models in family ([8]) are denoted by Gp^g, with < 
p < 2 and < q < 1, where p and q denote the maximum degree of t in and 
in a(t), respectively. The time-dependent GEV model ([7]) constructed in this way is 
a generalization of dH) (the latter is obtained by setting /ii = /i2 = cti = in ([8])). 
For model © with parameters ([8]), GEV inference amounts to estimating the parameter 
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vector 



(3 = [fio,fii,fi2,(To,ai,^] 



(9) 



by including time t as a covariate. 

Suppose we have a non- stationary dataset, from which a sequence of block maxima 
zt, with t = 1, . . . , m, is constructed. A log-likelihood function for the case ^ 7^ is 
defined as 



t=i 



/(/3) = -X] log + (1 + 1/0 log 



1+e 



a{t) 



1+e 



a{t) 



(10) 



(compare with provided that 



>0, z = l,...,m. 



(11) 



If ,^ = 0, an alternative lo g-likelihood function, derived from the Gumbel distribution, 



must be used (Coles 



20011). Numerical procedures are used to maximize the selected log- 



likelihood function, yielding the maximum likelihood estimate of the parameter vector 
/3. Con fidence inter vals for (3 may be computed by the expected or observed information 



matrix (Coles 



2001 ). 



Of course, all of the above procedure is performed in the spirit of "pure" inference. 
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that is determining the likelihood of the adopted parametric hypothesis and not its truth 
which, in the absence of a supporting theorem, remains unknown. Moreover, it should be 
kept in mind that several different models might fit the observations with similar reliabil- 
ity (likelihood). In this case, as no universal model is suggested or enforced (as opposed 
to the stationary case), there is no reason to prefer the one above the other. 



3c. Assessment of statistical models 



In the non-stationary context the analysis starts from a list of models (Gp ^ in our case, 
see dl])) which we fit to the data searching for the most adequate one. Assessment and 
comparisons of the inferences are based on standard graphical tools such as probabil- 
ity plot, quantile plot, and the likelihood ratio test. However, for the graphical model- 
checking the non-stationarity must b e taken into account. Reduction to Gumbel scale is 



2001 ) 



a practical way to treat this problem (|Colesl 

Let Zt,t= 1, . . . , m be a sequence of block maxima extracted from a non-stationary 
time series, from which the time-dependent GEV model G{'j2{t),a{t),^ ) has been fitted 
as described in the previous section. The sequence of maxima is transformed according 



to 



zt = log 



1 + e 



m. 



(12) 



If Zt are random variables with distribution G{'jl{t),a{t),^), then the transformation (fT2)) 
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produces variables Zt that have the standard Gumbel distribution: 



P{Zt < x) = exp (-e""^), X e M, 



(13) 



which is the GEV with parameters (/i, a, ^ = (0, 1, 0). Therefore, transformation (fT2)) 
attempts to remove the time-dependence from the sequence of maxima, bringing it as 
close as possible to the common scale given by the standard Gumbel distribution (fT3l) . 
This way, the distribution function and the quantiles of the transformed sequence of max- 
ima Zt can be compared with the empirical ones derived from (fT3l) . The probability plot 
is a graph of the pairs 



where Z(^i) < Z(^2) < ■ ■ ■ < Z(^rn) is the order statistics of the transformed sequence of 
maxima Zt. The quantile plot is given by the pairs 



For both plots, displacement of points from the diagonal indicates low quality of the 
inference. 

The likelihood ratio test is used to compare the goodness-of-fit of two nested models, 
that is, two models such that one of them is a sub-model (a particular case) of the other 




(14) 




(15) 
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one. Our family ^ of models is nested: for example Gi o is a sub-model of ^2,1^ 
obtained by setting 112 and ai to zero in the parameter vector /3 defined in Given 
integers < pi < p2 < 2 and < gi < g2 < 1, let /i and I2 be the maximized values of 
the log-likelihood (flOl ) for the nested models and Gp^^q^, respectively, and define 

the deviance statistic as 

V = 2{k-h}. (16) 

The likelihood ratio test states that the simpler model Gp^^g^ is to be rejected at the a- 
level of confidence in favor of Gp.^^^^ ifV > Ca, where Cq, is the (1 — a)-quantile of the 
xl distribution and k is the number of parameters that belong to Gp^^q^ and that are null 
in Gp^^g^ (k = 2 m our example above). 

Although the number of time-dependent models o ne may choo se from is virtually infi- 



nite, parsimony in the construction is reccommended (|Colesll2001l) : too many coefficients 
in the functions (/u(t), a{t)) would result in unacceptably large uncertainties, especially if 
few data are available. The search of the best model is carried out by trial-and-error: the 
choice of a more complex model should be strongly justified on theoretical grounds or by 
a significantly higher accuracy (that is, V exceedingly larger than Ca for nested models). 
However, the convergence of the above described procedure is by no means a guaran- 
tee of good estimate of the "true" probability density function: the latter is conceptually 
undefined. See our remarks at the end of Sec. [ 
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4. The time series: Total Energy of the Atmospheric Jet 



Model with a trend in average baroclinicity 



We consider here the same baroclinic jet model used in Part I, where the spectral or- 
der JT is set to 32. The model temperature is relaxed towards a given equator-to-pole 
profile which acts as baroclinic forcing. The statistical properties of the model radi- 
cally change when the parameter Te, determining the forced equator-to-pole tempera- 
tu re gradient, is varied. A phy si cal and dynamical de s cription of the model is given 



in 



Speranza and Malguzzil (|1988|) 



Malguzzi et al, 



(Il990h : 



Lucarini et al, 



(2006c 



In Part I we performed an extreme value analysis of the system's response with respect 
to variations in Te- Several stationary time series of the total energy E{t) were used as 
a basis for GEV inference. Each time series was generated with Te fixed at one value 
within a uniform grid on the interval [10, 50], with spacing of 2 units. We recall that, given 
the non-dimensionalization of the system, Te = I corresponds to 3.5 K, 1 unit of total 
energy corresponds to roughly 5 x 10^^ J, and t = 0.864 is one day, see Lucarini et al 



(|2006clld|) . In that case, all parameters of the system being kept fixed, after discarding an 



initial transient each time series of the total energy could be considered as a realization 
of a stationary stochastic process having weak long-range dependence. Therefore, the 
classical framework for GEV modeling was applied (see Sec.[3al). 

In the present setting, a specific linear trend is imposed on Te'- starting at time t = 0, 
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the model is run with a the time-dependent forcing parameter 

TE{t) = {T^-l)+tATE, te[0,to], (17) 

with = 10. Three values are chosen for the trend intensity ATe'- 2 units every Lb = 
1000, 300, and 100 years, yielding three time series for the total energy E(t). The range 
swept by TE{t) during integration is kept fixed in all three cases to the interval [9, 51], 
so that the total length of the time series depends on ATe- Each time series is split into 
21 data blocks B\ i = 1,...,21. The length Lb of each block corresponds to a time 
interval /* such that, as t varies within /*, the baroclinicity parameter TE{t) by (fTTI ) spans 
the interval 

[T^-1,T^ + 1], (18) 
which is 2 units wide and centered around one of the values considered in Part I: 

(TO , Ti, . . . , T|^) = (10, 12, 14, ... , 50). (19) 

Therefore, the total length L of the time series depends on the trend intensity, so that we 
have L = 21 X 2/ ATe = 21Lb. Moreover, since the time-span over which the maxima 
are computed is kept fixed to one year, the number of maxima in each data block B' also 
depends on ATe'. in fact, it is equal to Lb, see Tab.[T] 

Such a selection of the intervals as in (fTSi) allows for a direct comparison of the present 
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results with those obtained for stationary time series in Part 1. Moreover, our choices 
regarding block length and other factors are based on the indications provided in Part I, 
where the goodness-of-fit assessments performed by a variety of means showed that: 

• the adopted block length of one year makes sure that the extremes are uncorrelated 
and genuinely extreme; 

• the minumum length (100 data) used for the sequences of maxima yields robust 
inferences. 



5. Time-dependent GEV Analysis of the Total Energy 

For each data block B\ i = 1, . . . , 21, we first extract a sequence of yearly maxima 
zl, with t = 1, . . . , Lb- For compactness, each sequence is denoted in vector form as 
= (zl, zl, . . . , z\^). One GEV model of the form Gp^q (see ([8])) is fitted from each of 
the sequences 2'. For each i = 1, . . . , 21, the analysis follows three main steps: 

1. nested models Gp^q, for < p < 2 and < g < 1, are fitted on the i-th sequence 
of maxima z'; 

2. models with too low maximum likelihood are discarded and the deviance test is 
applied to the others to select the best estimate model; 

3. the best estimate model is graphically checked by examining the probability and 
quantile plots, and it is possibly rejected in favor of a model having less nonzero 
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parameters in the vector /3 as in (|9l). 

Following the above procedure, for each time interval I\ i = 1,...,21 time-dependent 
GEV models G^ qi{z'^) with parameters cr*(t), are inferred from the data block 

z\ Model Gpi qi{z^) (denoted for shortness G^ gi in the rest of this section) is the best 
estimate for the z-th data block, relative to the family of models Gp^g. Choosing a model 
with different orders (p, q) would either give poor results in the graphical checks, or fail 
to pass the likelihood ratio test. An example is given in Fig. [H for the data block i = 8 
in the time series with ATg = 2/(1000 years). The best estimate model has orders 
(p*, g*) = (1, 0). Models Gq^o and 6*2,1 are rejected, since they have too small likelihood 
and since the fit quality is very low, as it is illustrated by the probability and quantile 
plots. On the basis of the diagnostic plots, models Gi and Gi^i are both acceptable. 
However, the deviance statistic satisfies V = 2{/i_i — li^} = 3.64 < co.5 = 3.84 which 
is the 0.95-quantile of the xj-distribution. Therefore, as there is no strong support for 
selecting model Gi j, according to the likelihood ratio test the more parsimonious model 
Gi^o is preferred. 

Plots of the best estimate parameters (/i*(t), a*(t), ,^*) as functions of time are pro- 
posed in Fig. |2l Confidence intervals are computed as the worst case estimates: sup- 
pose that a model is chosen with p = 1, that is, for the best estimate /i(t) is linear 
/i*(t) = /iQ + jJL\t. Let a^i^ and cr^,. be the uncertainties in /Iq and respectively, pro- 
vided by the observed information matrix (see Sec. [S]). Then the confidence interval at 
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time t is computed as 

[p:\t) - 2(a^. + a^^^pit) + 2(a^. + a^^. (20) 

For most of the time intervals the best estimate model is such that and a^{t) are 
respectively linear and constant in time, that is, {p^, q^) = (1, 0). This is so for all models 
inferred with the fastest trend intensity ATe = 2/(100 years) (see Tab. IH), whereas for 
ATe = 2/(300 years) there are two exceptions: intervals i = 3 and i = 8, for which also 
a*(t) grows linearly in time {q' = 1, see Tab. [3]). For the slowest trend intensity we obtain 
= 1 for i = 2, 4, 5, 6, 7 and zero otherwise, whereas = 2 for i = 1, 2, 7 and = 1 
otherwise, see Tab.lH Summarizing, the best fits are mostly achieved by lowest order non- 
stationary models of the form {]?,q^) = (1,0). For slower trends, however, in some cases 
the best fit is of the form (p*, q^) = (1, 1) or even {j?, q^) = (2, 1). These cases typically 
occur for low i, that is, in the first portion of the time series. This is due to the fact that, 
for small Te (corresponding to small t through equation (fTTI)). although the hypothesis of 
smoothness, described in Sec. [21 may still be considered valid, the rate of variation of the 
SRB measure with respect to variations in Te is comparatively larger. To put it in simple 
words, the statistical behaviour (the attractor) of the baroclinic model is rather sensitive 
with respect to changes in Te. Therefore, the variation of the statistical properties in time 
is not quite "adiabatic", in the sense specified in Sec. [2l Correspondingly, a statistical 
model of enhanced complexity (more parameters) is needed to achieve goodness-of-fit. 
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In concluding this section, we emphasize that the convergence of the numerical pro- 
cedure used in the maximization of the likelihood function is here considerably more 
problematic than in the stationary case studied in Part 1. Indeed, in the present case it is 
often necessary to choose a good starting point for the maximization procedure in order 
to achieve convergence. For example, in several cases, after achieving convergence for a 
GEV model with order, say (p, q) = (2, 1), by using the inferred values of the parame- 
ter vector /3 in (|9l) as starting values for the maximization, a better fit (larger likelihood) 
of lower order is obtained. In fact, this has allowed us to reduce the total number of 
inferences with p = 2 and/or q = 1. 

6. Trend assessment 

When dealing with non-stationary data, the problem of assessing the sensitivity of trend 
inferences is particularly delicate. Beyond the serious conceptual problems explained in 
Sec. [21 one is confronted with several practical issues. Most of the sensitivity tests in 
Part I were based on examining a shorter portion of same time series or on calculating 
the maxima on data blocks of different lengths. In the present non- stationary context, 
both operations would result in an alteration of the statistical properties of the sample 
(exactly because of the non-stationarity) and this makes comparisons somewhat ambigu- 
ous. An example is provided in Fig. [31 where we compute the best estimate GEV fits 
using sequences of yearly maxima having different lengths - but starting at the same 
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instant (year 14500) - extracted from the time series with the slowest trend intensity 
ATe = 2/(1000 years). Notice that the best fit obtained by taking 100 yearly maxima 
is stationary. The corresponding extrapolations in time are, of course, completely wrong. 
By using 500 and 1000 maxima, the best estimates obtained (not shown) fall inside the 
confidence band of the 2000 years-based estimate for most values of time. 

The above example illustrates the trend dilemma: on the one hand, in order to be 
detected, a statistical trend has to be sufficiently fast with respect to the length of the 
record of observations; on the other hand, if the trend is too fast then the "adiabatic 
hypothesis" discussed in Sec.|2]is no longer valid: one is left with no "reference statistics" 
against which the inferred models can be compared. 

Moreover, when considering large time spans a further practical complication arises: 
due to the nonlinear dependence of the statistical properties with respect to the external 
parameter Te, a functional relation between the GEV parameters and time might require 
many parameters to achieve goodness-of-fit. Therefore, one faces the problem of large 
uncertainties in the parameter estimates or even lack of convergence. This has indeed 
been observed for the present time series: if we consider a long record, such that the 
change in Te is large, the model family Gp ^ with parameters as in ([8]) becomes inade- 
quate to catch the time dependence of the statistics of extremes. A first indication of this 
was reported at the end of Sec. [51 As a further example, we have examined a data block 
of length 5000 starting at year 14500 in the time series with ATg; = 2/(1000 years). In- 
spection of graphical diagnostics (probability and quantile plots) reveals that no model in 
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the family Gp ^ produces an acceptable inference. It should be emphasized that goodness- 
of-fit is achieved for the same time series using blocks of length 1000, that is, performing 
inferences that are more localized in time. So in this case the problem is not the failure of 
the "adiabatic" or smoothness hypothesis, but the nonlinear dependence of the attractor 
on the parameter Te, which manifests itself on sufficiently large time intervals. 



7. Smooth dependence on the forcing 

The set-up of the present analysis (see SecH]) has been chosen to allow comparison of the 
non-stationary GEV inferences with the results of Part I, obtained from statistically sta- 
tionary time series. To perform the comparison, for each i = 1, . . . , 21 the best estimate 
parameters (yu'(t), ,^*) inferred from data block are first expressed as functions of 
Te inside the interval (fTSi) . This is achieved by inverting the trend formula (flTl) . (writing 
time as a function of Te): 

t{TE) = '^''~l^^\ TeG [9,51]. (21) 

and inserting this into the expression of ^*). This yields functions that are 

denoted as (/i* (Tg) , a* (Tg) , . These are evaluated at the central point Tg of the interval 
of definition and plotted in Fig. IH Confidence intervals are given by the same worst 
case estimate (l20l) used for Fig. [21 A rather smooth dependence on Te is observed, 
especially for the GEV parameters fi and a. The location parameter fi turns out to be 
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not very sensitive to changes in the trend intensity, being much more sensibly dependent 
on variations in Te- Moreover its confidence intervals are always very small (relatively 
to the size of /x). 

Denote by 6*0,0 (i^*) the time-independent GEV models inferred in Part I from sta- 
tionary sequences of 1000 yearly maxima, computed with Te fixed at T^. Denote as 
/i\ a*, and the inferred values of the GEV parameters of GQfl{w^). Since the graphs 
of the parameters /i*, a\ and versus T^e ^^^Y closely match those in Fig.Hl comparison 
with the stationary data is presented under the form of relative differences (Fig. [5]). To be 
precise, on the left column the absolute values of the ratios 

are plotted against T|; (similarly for the GEV parameters cr and Remarkable agreement 
is obtained for the parameter /x: the relative differences less than 10% and drop below 5% 
for large Te and for all considered trend intensities. Excellent agreement is also obtained 
for cr (particularly for large Te) and for except for the fastest trend intensity ATe = 
2/(100 years). In the latter case, indeed, the sample uncertainty is as large as (or even 
larger than) the estimates self. 

We emphasize that inferring time-independent models G'o,o(^*) from the non-stationary 
data would induce very large errors (particularly in the scale and shape parameters), 
compare Fig. [B A much better (even surprising) agreement between the stationary and 
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non-stationary estimates is obtained with the procedure described in the previous section: 
first fitting the time-dependent model Gpi and then evaluating its parameters at the 
central point T^. There is agreement even in the estimates of the parameter which is 
usually the most difficult one to infer. Indeed, the inferred values are negative, in the case 
of stationary time series, since the attractor is bounded and since the energy observable 
E{t) is a continuous function of the phase space variables, the total energy is bounded on 
any orbit lying on (or converging to) the attractor. Therefore, the total energy extremes 
are necessarily WeibuU distributed is negative) Part I. Although this property is not 
bound to hold for non- stationary forcing, it is still verified, see Tab.[2l 

Two distinct power law regimes are identified for the GEV parameters (/i*, a*, ^*) as 
functions of T^, having the form 

fTiTh) = a^inr^ and a^Th) = a^i^y^ , (23) 

ee Fig. [6] and Fig. |7] The values of the fitted exponents 7^ and 70- in each scaling regime 
are reported in Tab. [5] and Tab. [6l respectively. A similar power law dependence of the 
GEV parameters on Te was already observed in Part I for the stationary data sets w': 
indeed, the exponents obtained there are very similar to those in Tab. |5] and Tab. [6l par- 
ticularly for large Te- The lack of a power-law scaling regime for the parameter a for 
small Te explains both the more pronounced differences between the stationary and non- 
stationary estimates (Fig.|5]) and the necessity of including a quadratic term in /i and/or a 
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linear term for o in the statistical model to get acceptable inferences. This highlights the 
strongly nonlinear behaviour of the baroclinic model, whose response to changes of Te 
has different features depending on the considered range of variation. 

Two factors explain the qualitative analogies and the quantitative agreements between 
the time-dependent models discussed here and the stationary results of Part I. First of all, 
the trend intensity imposed on in all cases is sufficiently slow with respect to the time 
of relaxation of the baroclinic model to the statistics of extreme values of the total en- 
ergy. For clarity, we emphasize that the latter time scale is that used in Sec. |2]to define 
the "adiabatic" hypothesis: it is the time necessary to obtain a good sampling of the SRB 
measure on the attractor, provided that one may consider the system as "frozen" (with 
constant Te) for sufficiently long time spans. We do not know whether this time scale 
bears any physical relation with other time scales, such a s those of baroclinic ins t ability 



or low-frequency variability (both have been described in lSperanza and Malguzzil (|1988|) 
for the present model). The second factor is that the system's statistical behaviour re- 
sponds rather smoothly to the imposed time-dependent variation of the parameter Te- 
This smooth depende nce on Te of t h e statis tical properties of the baroclinic model was 



analyzed in detail in iLucarini et al.l (|2006clldl) by considering not only global physical 



quantities such as total energy and average wind profiles, but also finer dynamical indi- 
cators, such as the Lyapunov exponents and dimension. Both properties of smoothness 
and "adiabaticity" are of crucial importance in order to justify the usage of non-stationary 
GEV models that are (locally) smooth functions of time, such as our polinomial family 
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8. Summary and Conclusions 

In this paper we have proposed a general, akhough not universal, framework for the anal- 
ysis of trends in extremes of climatic time series. When all the shortcomings which are 
present in datasets and observations have to be considered, a rigorous definition of ex- 
tremes and a neat, clean, and legible approach to the evaluation of trends is necessary 



in order to get useful and reliable information (IZhang et al.ll2005|) . The time-dependent 
approach allows to express the inferred GEV distributional parameters as functions of 
time. As expected, it is found that trend in the statistics of extreme values is detectable 
in a reliable way, provided that the record of observations is sufficiently long, depend- 
ing on the time scale of the trend itself. Trend inference and assessment is much more 
problematic than in the statistically stationary inference. First of all, one is faced with 
a serious conceptual problem: there is no "operational" definition of probability, since, 
to say it in loose words, the time series is not a sampling of a unique probability dis- 
tribution, as it is in the stationary case. Even if one assumes that the time series is a 
realization of a sequence of random variables (with different distributions), the statistical 
properties of the sample are altered by any operation such as resampling or taking shorter 
subsamples, which makes sensitivity studies somewhat ambiguous. One must assume 
that the distributions of the random variables vary slowly and smoothly with time, so that 
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the time series contains sufficient sampling information on the "local" (in time) statistical 
behaviour. 

In the present context, we have adopted GEV models whose parameters are poly- 
nomial in time: the location parameter /i is at most quadratic with respect to time and 
the scale parameter a is at most linear in time. Since the relation between the macro- 
scopic forcing Te and time is invertible, the time dependence of the inferred GEV mod- 
els can be expressed as a relation between the GEV parameters and Te, showing rather 
interesting properties. The location and scale parameters feature power-law dependence 
with respect to Te, while the shape parameter has in all cases a negative value. As ex- 
pected, both results are in agreement with what obtained in the companion paper (Part I) 
for stationary data. Since the parameter Te increases monotonically in the simulations 
with the baroclinic model, the system certainly does not possess any invariant measure. 
However, the results suggest that, as Te increases, the system explores statistical states 
which vary smoothly with Te and whose properties are locally quite similar to those ob- 
tained in the stationary setting. This is even captured for the relatively fast trend intensity 
/STe = 2/ (100 years). The proposed explanation is that: 



1 . the sys tem's statistical properties depend rather smoothly on Te (also compare (|Lucarini et al 



2006c 



MD); 

2. the adopted time-scales of variation of Te {i.e., the trend intensity AT^) are suffi- 
ciently slow compared to the relaxation time to the statistics of extreme values. 
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The second condition, that was explained in more detail in Sec. [2l amounts to the heuris- 
tic statement that for sufficiently short time spans the system's statistical properties can 
be considered /rozen to those holding for a corresponding value of Tg. The possibility 
of using GEV models that are locally smooth (polynomial functions of time) depends 
essentially on these two conditions. For example, for a system having several bifurca- 
tions as the control parameter is changed the time-dependent GEV modeling would be 
much more complicated. This problem is currently under investigation. However, even if 
the above two conditions do hold, the inference of time-dependent GEV models is valid 
locally in time, that is, if the sequences of maxima used for the inference span not too 
large time periods. For large time spans, indeed, the non-linear response of the baroclinic 
model to variations in Te becomes dominant and polynomial GEV models are no longer 
suitable. On the other hand, if the sequences of maxima used for the inference are too 
short (depending on the trend intensity), wrong trend estimates may be obtained. 

We conclude by observing that the present and the companion paper (Part I) are de- 
voted not merely to the statistical inference of extremes and their trends but also to explore 
the possibility of using extreme statistics in diagnosing the dynamical state of a geophys- 
ical fluid. Our analysis of the problem reveals, in fact, that diagnostics which is based 
on "universal" (GEV theorem), robust (smoothness properties), simple (power-law scal- 
ing), controllable (low-dimensional parametric) statistical mo dels can help very mu ch in 



setting up well targeted models of the general circulation, see (|Lucarini et al 



2006c 



There are several ways in which we plan to extend the present study. First of all, we 
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have considered a rather global indicator, the total energy of the system. Other choices 



might be to analyze the wave kinetic energy, the available energy or also the maximum 
vorticity on the domain of the model, which might behave differently as Te is changed. 
Moreover, there are delicate issues connected with reducing the scale from a global indi- 
cator to a local one, such as the value of the wind on a grid point. This brings into play 
all complications due to the multifractality and the spatial dependence of the process. A 
further development of the present work is the us age of extreme statistics as a dynamical 



indicator, in the sense of process oriented metrics 



Lucarini et al 



(l2006a|) . All these issues 



are currently under investigation. 
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List of Figures 

1 Diagnostic plots of GEV inferences with model ^ and parameters as 
in ([8]), for block Bi with i = 8 (corresponding to = 24) and ATe = 
2/(1000 years). Top row: probability plots. Bottom row: quantile plots. 
From left to right column: plots for models Gp^q (see ([8])), with (p, q) and 

the corresponding log-likelihood / (see Q) indicated on top |41 

2 GEV parameters as functions of time, for the three considered values 
of trend intensity ATe: from top to bottom, ATe = 2/(1000 years), 
2/(300 years), and 2/(100 years), respectively. For each trend intensity 
the inferred time-dependent parameters (/i(t), a{t),^) (left, center, right 
column, respectively) of the best estimate model Gpi^gi (2*) are plotted. . |42 

3 Parameter /i(t) of the best estimate GEV inferences as a function of time. 
The time series with slowest trend intensity ATe = 2/(1000 years) has 
been used, taking yearly maxima over a data block starting at year 14001. 
For legibility, only the confidence intervals have been plotted. Left: in- 
ferences obtained with 100, 500, and 1000 yearly maxima. The best 
estimate fit based on 100 data is stationary (/ii = 0) and it has been 
extrapolated to 500 years. Right: inferences obtained with 500, 1000, 
and 2000 yearly maxima. In all cases, the best estimate GEV model has 

p = 1, that is, n{t) = Ho + l^it is a linear function of time |43 
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Parameters {'fj^{T^),a^{T^),^^) (from left to right, respectively) of the 
best estimate GEV model ^^(z*) evaluated at the central point 
of each of the 21 intervals (fTS]) . From top to bottom the trend inten- 
sity ATe is equal to 2/(1000 years), 2/(300 years), and 2/(100 years), 

respectively 

Same as Fig. |4] for the estimates obtained with the stationary data in Part 
I, see text for details. The time-dependent estimates of Fig. |4] are plotted 

with dotted lines 

Power law fits of the inferred values of /?(Tg) as a function of (see (fT9l )) 
From left to right: trend intensities of 2/(1000 years), 2/(300 years), and 
2/(100 years) have been used. In each case, there are two intervals of Te 
characterized by different scaling law, separated by a point T|., compare 

Tab. El 

Same as Fig.[6]for the inferred values a^{T'^) 
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Figure 1 : Diagnostic plots of GEV inferences with model (|7]) and parameters as in ([8]), 
for block Bi with i = 8 (corresponding to = 24) and ATe = 2/(1000 years). Top 
row: probability plots. Bottom row: quantile plots. From left to right column: plots 
for models Gp^q (see ([8])), with (p, q) and the corresponding log-likelihood / (see 
indicated on top. 
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Figure 2: GEV parameters as functions of time, for the three considered values of 
trend intensity ATg;: from top to bottom, AT^ = 2/(1000 years), 2/(300 years), and 
2/(100 years), respectively. For each trend intensity the inferred time-dependent param- 
eters {fi{t),a{t),^) (left, center, right column, respectively) of the best estimate model 
G^^^{z^) are plotted. 
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Figure 3: Parameter fi{t) of the best estimate GEV inferences as a function of time. 
The time series with slowest trend intensity ATe = 2/(1000 years) has been used, taking 
yearly maxima over a data block starting at year 14001 . For legibility, only the confidence 
intervals have been plotted. Left: inferences obtained with 100, 500, and 1000 yearly 
maxima. The best estimate fit based on 100 data is stationary (/ii = 0) and it has been 
extrapolated to 500 years. Right: inferences obtained with 500, 1000, and 2000 yearly 
maxima. In all cases, the best estimate GEV model has p = 1, that is, fi{t) = + fxit is 
a linear function of time. 
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Figure 4: Parameters (ji\T^),a\T^), (from left to right, respectively) of the best 
estimate GEV model Gp^ q^{z^) evaluated at the central point T'^ of each of the 21 in- 
tervals (fTSi) . From top to bottom the trend intensity ATe is equal to 2/(1000 years), 
2/(300 years), and 2/(100 years), respectively. 
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Figure 5: Same as Fig.|4]for the estimates obtained with the stationary data in Part I, see 
text for details. The time-dependent estimates of Fig. |4] are plotted with dotted lines. 
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Figure 6: Power law fits of the inferred values of 'j2\T^) as a function of (see (fT9l)). 
From left to right: trend intensities of 2/(1000 years), 2/(300 years), and 2/(100 years) 
have been used. In each case, there are two intervals of Te characterized by different 
scaling law, separated by a point T^, compare Tab.|5l 
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Figure 7: Same as Fig.|6]for the inferred values a''{T^). 
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ATg 2/1000 2/300 2/100 
L 21000 6300 2100 
Lb 1000 300 100 

Table 1 : The length L of each of the three time series and the length Lb of each of the 21 

the data blocks Bi (both are expressed in years), as a function of the intensity ATg of the 

trend (flTl) imposed on the parameter Te of the baroclinic model. 
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Table 2: Best estimate GEV fits Gpi^qi{z^) with parameter vector as in Q for the non- 
stationary time series with trend intensity AT^ = 2/(1000 years), see text for details. 
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Table 3: As in Tab. |2] for trend intensity ATe = 2/(300 years). 
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Table 4: As in Tab.|2]for trend intensity ATe = 2/(100 years). 
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Table 5: Power law fits of the inferred location parameter /i*(T|;) as a function of 
(see (fT9l) ) of the form /i*(T|) = a^(T|)'''f'. Two distinct scaling regimes (with distinct 
exponents 7^ i and 7^ 2) are identified and the corresponding adjacent intervals in the 
T^-axis are separated by Tl. 
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Table 6: Same as Tab.[5]for the inferred scale parameter a\T^). 
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Abstract 

A baroclinic model for the atmospheric jet at middle-latitudes is used 
as a stochastic generator of non-stationary time series of the total energy of 
the system. A linear time trend is imposed on the parameter Te, descriptive 
of the forced equator-to-pole temperature gradient and responsible for set- 
ting the average baroclinicity in the model. The focus lies on establishing 
a theoretically sound framework for the detection and assessment of trend 
at extreme values of the generated time series. This problem is dealt with 
by fitting time-dependent Generalized Extreme Value (GEV) models to se- 
quences of yearly maxima of the total energy. A family of GEV models is 
used in which the location jj, and scale parameters a depend quadratically 
and linearly on time, respectively, while the shape parameter is kept con- 
stant. From this family, a model is selected by using diagnostic graphical 
tools, such as probability and quantile plots, and by means of the likelihood 
ratio test. The inferred location and scale parameters are found to depend in 
a rather smooth way on time and, therefore, on Te- In particular, power-law 
dependences of fi and a on Te are obtained, in analogy with the results of 
a previous work where the same baroclinic model was run with fixed values 
of Te spanning the same range as in this case. It is emphasized under which 
conditions the adopted approach is valid. 

PACS: 02.50.Tt, 02.70.-c, 47.1 l.-j, 92.60.Bh, 92.70.Gt 
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1. Introduction 

In the context of Climate Change, an intensely debated question is whether the statistics 
of extreme meteo-climatic events is changing (and/or will change) and, in case, how fast 
it is changing (and/or will change). For example, the role of time-dependence in the 
statistics of extreme weather events has been at the heart of discussions about climate 
change since the work by ?. In particular, the detection of trends in the frequency of 
intense precipitation has been the object of much research, particularly at regional level, 
see e.g. ?? for the USA and ?? for the Mediterranean area. The general relevance of the 
problem has been highlighted in the 2002 release of a specific IPCC report on Changes in 
extreme weather and climate events (available at http://www.ipcc.ch/pub/support.htm). In 
fact, the emphasis laid on the subject by the IPCC report reverberated in many countries 
the question "is the probability of major impact weather increasing?". This question 
reached the big public almost everywhere and innumerable studies of trends in series of 
"extremes" were undertaken. These studies mainly deal with variables of local character, 
typically precipitation and temperature at specific stations. Moreover, most studies are 
regional: see e.g. the proceedings of the Italia-USA meeting held in Bologna in 2004 (?) 
for the relevance of the extreme events in the Mediterranean Climates and the INTERREG 
IIIB - CADSES project HYDROCARE, http://www.hydrocare-cadses.net, 
for impacts of extreme events in the hydrological cycle of the central-eastern Europe. 
In a preceding, companion paper (?) (which we refer to by Part I in the sequel) 
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we have addressed the problem of extreme value statistical inference on statistically sta- 
tionary time series produced by a dynamical system providing a minimal model for the 
dynamics of the mid-latitudes baroclinic jet. There reported is, from mathematical lit- 
erature, a suitable, rigorous, "universal" setting for the analysis of the extreme events in 
stationary time series. This is based on Gnedenko's theorem (?) according to which the 
distribution of the block-maxima of a sample of independent identically distributed vari- 
ables converges, under fairly mild assumptions, to a member of a three-real parameter 
family of distributions, the so-called Generalized Extreme Value (GEV) distribution (?). 
The GEV approach to the analysis of extremse requires that three basic conditions are 
met, namely the independence of the selected extreme values, the consideration of a sujfi- 
ciently large number of extremes, the selection of values that are genuinely extreme. This 
could be performed relatively easily for the case at hand. 

Part I was originally motivated by the interest in weather having "major impact" (on 
human life and property) in the Mediterranean area, in particular intense precipitation and 
heat waves over Italy. See, for example, ??????? and the MEDEX Phase 1 report (avail- 
able at http : / /medex . inm . uib . es /) for related results and activities. The study 
reported in Part I has revealed, among other things, that diagnostics of extreme statis- 
tics can highlight interesting dynamical properties of the analyzed system. Properties 
which, thanks to the "universality of the GEV", can be investigated in a low dimensional- 
ity space of parametric probability density functions, although at the expenses of the total 
length of the observational record of the system in order to capture a sufficient number of 
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independent extremes. A key role (that is presently being explored elsewhere, in the con- 
text of general atmospheric circulation theory) was played in Part I by the smoothness 
of variation of the extreme statistics parameters (average, variance, shape factor) upon 
the external (forcing) parameters of the system. In this paper, again, we devote attention 
to exploring the statistics of extremes as a dynamical indicator, this time in the frame- 
work of the (typically meteorological) statistical inference problem of detecting trends in 
observations. 

The definition of a rigorous approach to the study of extremes is much harder when the 
property of stationarity does not hold. One basic reason is that there exists no universal 
theory of extreme values (such as e.g. a generalization of Gnedenko's theorem) for non- 
stationary stochastic processes. Moreover, in the analysis of observed or synthetically 
generated sequences of data of finite length, practical issues, such as the possibility of 
unambiguously choosing the time scales which defines the statistical properties and their 
changes, become of critical importance. Nevertheless, GEV-based statistical modeling 
offers a practical unified framework also for the study of extremes in non- stationary time 
series. In the applications, the three parameters of the GEV distribution are taken as time- 
dependent and time is introduced as a covariate in the statistical inference procedure (?). 
The practical meaning of this assumption is that the probability of occurrence (chance) of 
the considered extreme events evolves in time pretty much as we are inclined to think in 
our everyday life. However, giving a scientific meaning to such an assumption is possible 
only in an intuitive, heuristical fashion: in an "adiabatic" limit of infinitely slow trends 
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(but rigorously not even in such a limit). We adopt this point of view not only because 
it is in line with the common practice and view of extremes, but also because interesting 
dynamical properties can be inferred from extremes, in analogy with the findings in Part 
I. 

In the present paper we perform and assess time-dependent GEV inference on non- 
stationary time series E(t) of the total energy obtained by the same simplified quasi- 
geostrophic model that was used in Part I. The model undergoes baroclinic forcing to- 
wards a given latitudinal temperature profile controlled by the forced equator-to-pole 
temperature difference T^; see ?? for a thorough description. We analyze how the pa- 
rameters of the GEV change with time when a linear trend is imposed on the large scale 
macroscopic forcing Te, that is, when Te is taken as a (linear) function of time. Since 
this functional relation is invertible, we derive a parametrization relating the changes in 
the GEV to variations in Te (instead of time). One major goal here is to present a method- 
ological framework to be adopted with more complex models and with data coming from 
observations, as well as an assessment of the performance of the GEV approach for the 
analysis of trends in extremes in the somewhat grey area of non-stationary time series. 
Methodologically, our set-up is somewhat similar to that of ? regarding the procedures 
of statistical inference. However, in this case we face two additional problems: 

1. as in Part I, the statistical properties of the time series E{t) cannot be selected 
a priori: in the stationary case {Te constant in time) and much less in the non- 
stationary case there is no explicit formula for the probability distribution of the 
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observable E(t); 

2. moreover, in the non-stationary case we even lack a definition (and in fact a mere 
candidate) of what might be the probability distribution of E(t): certainly not a 
frequency limit for t ^ cxd (and not by construction, as opposed to ? who use 
genuinely stochastic generators). 

This also means that we have no hypothesis concerning the functional form of the trend in 
the statistics of extremes of E(t), resulting from the trend imposed on the control param- 
eter Te. The lack of a general GEV theorem for non- stationary sequences implies that 
the choice of the time-dependent GEV as a statistical model is, in principle, arbitrary: 
other models might be equally (or better) suitable. Here comes into play the "adiabatic- 
ity" hypothesis mentioned above, which leads to the central statement of this paper: if 
the trend is sufficiently slow and if the statistical behaviour of the atmospheric model has 
a sufficiently regular response with respect to variations in the external parameters, the 
GEV remains a suitable model for inference of trend in extremes. 

The structure of the paper follows. In Sec. ?? we describe the general problem of 
the characterization of statistical trends in deterministic models, with both its conceptual 
and practical implications. Then in Sec. ?? we describe how the GEV modeling can be 
applied to non- stationary time-series and how the quality-check of the fits is performed. 
In section Sec. ?? we present the time series considered in this work and the set-up of 
the numerical experiments performed with the atmospheric model. The inferences for 
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various values of the trend in the forcing parameter Te are presented in Sec. ?? and a 
sensitivity analysis is carried out in Sec. ??. Comparison with the stationary case ana- 
lyzed in Part I is given in Sec. ??. In Sec. ??, finally, we summarize the main findings of 
this work, highlighting the future research developments. 



The stochastic generator used in this paper to produce the time series is a deterministic 
model (an ordinary differential equation), whose dynamics, for the considered range of 
values of Tg, is chaotic in the sense that it takes place on a strange attractor A in phase 
space (?). See ?? for a study of the properties of this attractor, including sensitivity 
with respect to initial conditions. The statistical behaviour of this type of time series is 
determined by the Sinai-Ruelle-Bowen (SRB) probability measure /i (?): this is a Borel 
probability measure in phase space which is invariant under the flow /* of the differential 
equation, is ergodic, is singular with respect to the Lebesgue measure in phase space 
and its conditional measures along unstable manifolds are absolutely continuous, see ? 
and references therein. Moreover, the SRB measure is also a physical measure: there 
is a set V having full Lebesgue measure in a neighbourhood t/ of A such that for every 
continuous observable (j) : U ^ M, we have, for every x ^ V, the frequency-limit 
characterization 



2. Statistical trends: the theoretical problem 




(1) 
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Existence and uniqueness of an SRB measure ji have been proved only for very special 
classes of flows /* (in particular, for flows that possess an Axiom-A attractor, see ?). 
However, existence and uniqueness of /x are necessary to define a stationary stochastic 
process associated to an observable (f). In turn, this allows to consider a given time series 
of the form cjy^f^^x)) : t > as a realization of the stationary process, justifying statistical 
inference on a solid theoretical basis. In part I, we conjectured existence an uniqueness 
of an SRB measure for the atmospheric model, providing the theoretical foundation to 
the application of GEV models in the inference of extreme values. 

In certain cases of non-autonomous ordinary differential equations (for example, if the 
dependence on time is periodic), it still is possible to define, at least conceptually, what 
an SRB measure is. However, in the present case, due to the form of time-dependence 
adopted for the parameter Te, the atmospheric model admits no invariant measure. This 
means that there is no (known) way to associate a stochastic process to the time series 
of the total energy. In other words, it is even in doubt what we mean by "statistical 
properties" of the time series, since it is impossible to define a probability distribution. 
This conceptual problem has a very serious practical consequence: the "operational" 
definition of probability as a frequency limit (as in (??)) is not valid in the non- stationary 
case, since the time evolution is not a sampling of a unique probability distribution. Even 
if one assumes the existence of a sequence of distinct probability distributions, one for 
each observation, one realization (the time series) does not contain sufficient statistical 
information, since each distribution is very undersampled (with only one observation). 
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Despite all these problems, the results in Part I suggest a framework which is, for 
the moment, formulated in a heuristic way. Suppose you evolve an initial condition x in 
phase space by the flow /* of the autonomous atmospheric model, that is the system in 
which Te is kept fixed to some value T^. After an initial time span (transient), say for t 
larger than some to > 0, the evolution may be thougth to take place on the attractor 
A and time-averages of the form 



which is the average of by the SRB measure /iq existing at the value Te = (the 
"attractor average" at T^). Now suppose that at some ti > to the parameter Te is abruptly 
changed to some value > T^: there will be some transient interval, call it [ti,t2], 
during which the evolution approaches the new attractor of the atmospheric model, 
that is the attractor existing for Te fixed at T^. After that time span, the evolution may 
be considered to take place on the new attractor. 

In our case, though Te varies continuously (linearly) with time, if the trend magnitude 
is low, then Te may be considered constant (with good approximation) during time spans 




(2) 



may be considered as approximations of 




(3) 
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that are sufficiently long in order to have both convergence to the "new" attractor and good 
sampling of the "new" SRB measure, in the sense sketched above. Though admittedly 
heuristic, this scenario allows to clarify under which condition it is still reasonable to 
speak of "statistical properties" of a time series generated by a non-autonomous system: 
namely, the closeness to a stationary situation for time spans that are sufficiently long. 
This is the "adiabatic" hypothesis which we mentioned in the introduction. An essential 
ingredient for this to hold is that the statistical properties of the autonomous model do 
not sensibly depend on the external parameter Te, in the sense that no abrupt transitions 
(bifurcations) should take place as Te is varied. This was indeed checked for the system 
at hand in Part I. Notice that the validity of the "adiabatic" hypothesis also has a useful 
practical consequence: one can use the statistics of the stationary system as a reference 
against which the results for the non-stationary case can be assessed. Having this scenario 
in mind, we proceed to the description of the time-dependent GEV approach in the next 
section. 



3. GEV modelling for non-stationary time series 

3a. Stationary case 

The GEV approach for sequences of independent, identically distributed (i.i.d.) random 
variables is by now rather standard (??????????). The foundation is Gnedenko's theo- 
rem (??). Consider a time series, assumed to be a realization of a sequence of i.i.d. ran- 
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dom variables. The time series is divided into m consecutive time-frames (data blocks), 
each containing n subsequent observations, equally spaced in time. Denote hy zi, . . . , Zm 
the sequence of the block maxima taken over each time-frame. Under fairly mild as- 
sumptions, the distribution of the block-maxima converges, in a suitable limit involving 
a rescaling, to the GEV distribution, defined as 



G{x; fi, a, = exp 



1+e 



X — jJ 



a 



(4) 



for all X in the set {x : 1 + S,{x — ^l)/o > 0} and G{x) = otherwise, where a > and 
^ 7^ 0. If ^ = the limit distribution is the Gumbel distribution 



G{x; /i, cr, 0) = exp I — exp 



X — fl 



a 



X G 



(5) 



GEV inference consists in estimating the distributional parameters (/i, cr, ^) (called lo- 
cation, scale and shape parameter, respectively) from the available data. A widely used 
technique consists in numerically maximizing the log-likelihood function a, For 
^ 7^ 0, this is defined as 



— m log a — \ \ 



1+^ 



2t - /i 



cr 



1+e 



Zt- 



a 



(6) 
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while an analogous formula holds for ^ = (?). 

A generalization of the GEV theorem holds for time series that are realizations of 
stationary stochastic processes such that the long-range dependence is weak at extreme 
levels (??). In the applications, this property is assumed to hold whenever the block 
maxima are uncorrected for sufficiently large block sizes. In this case, GEV inference 
and assessment is carried out by the same tools (maximum likelihood, diagnostic plots, 
etc.) used in the i.i.d. context, see ?. 

3b. Time-dependent case 

If stationarity of the time series does not hold, then the limiting distribution function is 
no longer bound to be the GEV or any other family (??): no theories of extreme values 
exist in this context. Some exact results are known only in certain very specialized types 
of non- stationarity (??), but it is very unlikely that a general theory can be established. 
However, GEV-based statistical modeling of extreme values can be performed also in 
the case of time-dependent phenomena by adopting a pragmatic approach, where the 
GEV distribution (??) is used as a template: time-dependent parameters and cr(t) 
are considered, yielding a GEV model of the form 

G{x■^^{t),a{t),i). (7) 
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Usually ^ is kept time-independent in order to avoid numerical problems, since it is usu- 
ally the most delicate parameter to estimate (?). Different kinds of time-dependence can 
be imposed for {n{t), (j{t)). In this paper, we adopt a simple polynomial family of mod- 
els: 



with /io,i,2 and ctq,! € M. GEV models in family (??) are denoted by Gp^g, with < 
p < 2 and < g < 1, where p and q denote the maximum degree of t in and 
in a{t), respectively. The time-dependent GEV model (??) constructed in this way is a 
generalization of (??) (the latter is obtained by setting /ii = /i2 = ai = in (??)). For 
model (??) with parameters (??), GEV inference amounts to estimating the parameter 
vector 



by including time t as a covariate. 

Suppose we have a non- stationary dataset, from which a sequence of block maxima 
Zt, with t = 1, . . . , m, is constructed. A log-likelihood function for the case ,^ 7^ is 
defined as 



H{t) = /io + filt + /i2t^, 



a{t) = (To + (Tit, 



(8) 



/3 = [/^o,/^i,/^2,o-o,o-i,^] 



(9) 



m f 



t=i ^ 



( 



Zt - jJ'jt) 

a{t) 



) 



+ 




(10) 
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(compare with (??)), provided that 

>0, z = l,...,m. (11) 

If ,^ = 0, an alternative log-likelihood function, derived from the Gumbel distribution, 
must be used (?). Numerical procedures are used to maximize the selected log-likelihood 
function, yielding the maximum likelihood estimate of the parameter vector (3. Confi- 
dence intervals for /3 may be computed by the expected or observed information ma- 
trix (?). 

Of course, all of the above procedure is performed in the spirit of "pure" inference, 
that is determining the likelihood of the adopted parametric hypothesis and not its truth 
which, in the absence of a supporting theorem, remains unknown. Moreover, it should be 
kept in mind that several different models might fit the observations with similar reliabil- 
ity (likelihood). In this case, as no universal model is suggested or enforced (as opposed 
to the stationary case), there is no reason to prefer the one above the other. 

3c. Assessment of statistical models 

In the non-stationary context the analysis starts from a list of models (Gp ^ in our case, 
see (??)) which we fit to the data searching for the most adequate one. Assessment and 
comparisons of the inferences are based on standard graphical tools such as probabil- 
ity plot, quantile plot, and the likelihood ratio test. However, for the graphical model- 
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checking the non-stationarity must be taken into account. Reduction to Gumbel scale is 
a practical way to treat this problem (?). 

Let zt,t = 1, . . . , m be a sequence of block maxima extracted from a non-stationary 
time series, from which the time-dependent GEV model G{'j2{t),a{t),^ ) has been fitted 
as described in the previous section. The sequence of maxima is transformed according 
to 

t = l,...,m. (12) 

If Zt are random variables with distribution G{'jl{t),a{t),$^), then the transformation (??) 
produces variables Zt that have the standard Gumbel distribution: 

P{Zt <x)= exp (-e-^), X G M, (13) 



Zt = log 



a{t) 



which is the GEV with parameters (/i, a, ^) = (0, 1, 0). Therefore, transformation (??) 
attempts to remove the time-dependence from the sequence of maxima, bringing it as 
close as possible to the common scale given by the standard Gumbel distribution (??). 
This way, the distribution function and the quantiles of the transformed sequence of max- 
ima it can be compared with the empirical ones derived from (??). The probability plot 
is a graph of the pairs 



^ ;exp(-e-^(^)) 1 , j = l,...,m, (14) 



m + 1 
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where 2(1) < zi^2) < ■ • ■ < is the order statistics of the transformed sequence of 
maxima The quantile plot is given by the pairs 

(-log(-log(;;7^)) i ' J = l,...,m. (15) 

For both plots, displacement of points from the diagonal indicates low quality of the 
inference. 

The likelihood ratio test is used to compare the goodness-of-fit of two nested models, 
that is, two models such that one of them is a sub-model (a particular case) of the other 
one. Our family ^ of models is nested: for example Gi is a sub-model of ^2,1, 
obtained by setting fi2 and ai to zero in the parameter vector /3 defined in (??). Given 
integers < < p2 < 2 and < gi < g2 < let /i and I2 be the maximized values of 
the log-likelihood (??) for the nested models Gp^^g^ and Gp^^q^, respectively, and define 
the deviance statistic as 

V = 2{k~h}. (16) 

The likelihood ratio test states that the simpler model Gp^^gj is to be rejected at the a- 
level of confidence in favor of Gp2 ^3 ifV > Ca, where Cq, is the (1 — a) -quantile of the 
xl distribution and k is the number of parameters that belong to Gp2,q2 ^iid that are null 
in Gp^^gj (k = 2'm our example above). 

Although the number of time-dependent models one may choose from is virtually 
infinite, parsimony in the construction is reccommended (?): too many coefficients in 
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the functions cr(t)) would result in unacceptably large uncertainties, especially if 

few data are available. The search of the best model is carried out by trial-and-error: the 
choice of a more complex model should be strongly justified on theoretical grounds or by 
a significantly higher accuracy (that is, V exceedingly larger than Cq, for nested models). 
However, the convergence of the above described procedure is by no means a guarantee 
of good estimate of the "true" probability density function: the latter is conceptually 
undefined. See our remarks at the end of Sec. ??. 



4. The time series: Total Energy of the Atmospheric Jet 
Model with a trend in average baroclinicity 

We consider here the same baroclinic jet model used in Part I, where the spectral order 
JT is set to 32. The model temperature is relaxed towards a given equator-to-pole profile 
which acts as baroclinic forcing. The statistical properties of the model radically change 
when the parameter Te, determining the forced equator-to-pole temperature gradient, is 
varied. A physical and dynamical description of the model is given in ????. 

In Part I we performed an extreme value analysis of the system's response with respect 
to variations in Te- Several stationary time series of the total energy E(t) were used as 
a basis for GEV inference. Each time series was generated with Te fixed at one value 
within a uniform grid on the interval [10, 50], with spacing of 2 units. We recall that, given 
the non-dimensionalization of the system, = 1 corresponds to 3.5 K, 1 unit of total 
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energy corresponds to roughly 5 x 10 J, and t = 0.864 is one day, see ??. In that case, all 
parameters of the system being kept fixed, after discarding an initial transient each time 
series of the total energy could be considered as a realization of a stationary stochastic 
process having weak long-range dependence. Therefore, the classical framework for 
GEV modeling was applied (see Sec. ??). 

In the present setting, a specific linear trend is imposed on T^: starting at time t = 0, 
the model is run with a the time-dependent forcing parameter 

TE{t) = {T^-l)+tATE, te[0,to], (17) 

with = 10. Three values are chosen for the trend intensity AT^: 2 units every Lb = 
1000, 300, and 100 years, yielding three time series for the total energy E{t). The range 
swept by Tj;(t) during integration is kept fixed in all three cases to the interval [9, 51], 
so that the total length of the time series depends on ATe- Each time series is split into 
21 data blocks B\ i = 1, . . . , 21. The length Lb of each block corresponds to a time 
interval P such that, as t varies within /*, the baroclinicity parameter TE{t) by (??) spans 
the interval 

[Ti;-1,T^ + 1], (18) 
which is 2 units wide and centered around one of the values considered in Part I: 

(T°,ri,...,T#) = (10,12,14,...,50). (19) 
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Therefore, the total length L of the time series depends on the trend intensity, so that we 
have L = 21 X 2/ATe = 2\Lb- Moreover, since the time- span over which the maxima 
are computed is kept fixed to one year, the number of maxima in each data block 5* also 
depends on AT^;: in fact, it is equal to Lb, see Tab. ??. 

Such a selection of the intervals as in (??) allows for a direct comparison of the present 
results with those obtained for stationary time series in Part I. Moreover, our choices 
regarding block length and other factors are based on the indications provided in Part I, 
where the goodness-of-fit assessments performed by a variety of means showed that: 

• the adopted block length of one year makes sure that the extremes are uncorrected 
and genuinely extreme; 

• the minumum length (100 data) used for the sequences of maxima yields robust 
inferences. 

5. Time-dependent GEV Analysis of the Total Energy 

For each data block _B*, i = 1, . . . , 21, we first extract a sequence of yearly maxima 
zl, with t = 1, . . . , Lb- For compactness, each sequence is denoted in vector form as 
z'^ = {z{, ■ ■ . , z\^). One GEV model of the form Gp^q (see (??)) is fitted from each of 
the sequences 2'. For each i = 1, . . . , 21, the analysis follows three main steps: 

1. nested models Gp^q, for < p < 2 and < g < 1, are fitted on the i-th sequence 
of maxima 2 '; 
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2. models with too low maximum likelihood are discarded and the deviance test is 
applied to the others to select the best estimate model; 

3. the best estimate model is graphically checked by examining the probability and 
quantile plots, and it is possibly rejected in favor of a model having less nonzero 
parameters in the vector /3 as in (??). 

Following the above procedure, for each time interval P, i = 1,...,21 time-dependent 
GEV models G^ff{z'') with parameters cr*(t), are inferred from the data block 

2*. Model Gpi^gi{z^) (denoted for shortness 5. in the rest of this section) is the best 
estimate for the i-th data block, relative to the family of models Gp^g. Choosing a model 
with different orders {p, q) would either give poor results in the graphical checks, or 
fail to pass the likelihood ratio test. An example is given in Fig. ??, for the data block 
z = 8 in the time series with ATg = 2/(1000 years). The best estimate model has orders 
(p*, q") = (1, 0). Models Go,o and G2,i are rejected, since they have too small likelihood 
and since the fit quality is very low, as it is illustrated by the probability and quantile 
plots. On the basis of the diagnostic plots, models Gi^ and Gi^i are both acceptable. 
However, the deviance statistic satisfies V = 2{/i^i — /i_o} = 3.64 < co.5 = 3.84 which 
is the 0.95-quantile of the Xi-distribution. Therefore, as there is no strong support for 
selecting model Gi 1, according to the likelihood ratio test the more parsimonious model 
Gi^o is preferred. 
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Plots of the best estimate parameters as functions of time are pro- 

posed in Fig. ??. Confidence intervals are computed as the worst case estimates: sup- 
pose that a model is chosen with p = 1, that is, for the best estimate ii{t) is linear 
= /ig + ii\t. Let a^i^ and cr^i be the uncertainties in /Iq and respectively, pro- 
vided by the observed information matrix (see Sec. ??). Then the confidence interval at 
time t is computed as 

- 2(a^. + v^t), + 2(a^. + a^^t)]. (20) 

For most of the time intervals the best estimate model is such that and are 
respectively linear and constant in time, that is, {j?, q^) = (1, 0). This is so for all models 
inferred with the fastest trend intensity ATe = 2/(100 years) (see Tab. ??), whereas for 
ATe = 2/(300 years) there are two exceptions: intervals i = 3 and i = 8, for which 
also a*(t) grows linearly in time (g* = 1, see Tab. ??). For the slowest trend intensity we 
obtain = 1 for i = 2, 4, 5, 6, 7 and zero otherwise, whereas = 2 for 2 = 1,2,7 and 
= 1 otherwise, see Tab. ??. Summarizing, the best fits are mostly achieved by lowest 
order non-stationary models of the form (p*, q^) = (1, 0). For slower trends, however, in 
some cases the best fit is of the form (p^, q^) = (1, 1) or even (p*, q^) = (2, 1). These 
cases typically occur for low i, that is, in the first portion of the time series. This is due 
to the fact that, for small Te (corresponding to small t through equation (??)), although 
the hypothesis of smoothness, described in Sec. ??, may still be considered valid, the 
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rate of variation of the SRB measure with respect to variations in Te is comparatively 
larger. To put it in simple words, the statistical behaviour (the attractor) of the baroclinic 
model is rather sensitive with respect to changes in Te. Therefore, the variation of the 
statistical properties in time is not quite "adiabatic", in the sense specified in Sec. ??. 
Correspondingly, a statistical model of enhanced complexity (more parameters) is needed 
to achieve goodness-of-fit. 

In concluding this section, we emphasize that the convergence of the numerical pro- 
cedure used in the maximization of the likelihood function is here considerably more 
problematic than in the stationary case studied in Part I. Indeed, in the present case it is 
often necessary to choose a good starting point for the maximization procedure in order 
to achieve convergence. For example, in several cases, after achieving convergence for a 
GEV model with order, say (p, g) = (2, 1), by using the inferred values of the parameter 
vector /3 in (??) as starting values for the maximization, a better fit (larger likelihood) 
of lower order is obtained. In fact, this has allowed us to reduce the total number of 
inferences with p = 2 and/or q = I. 

6. Trend assessment 

When dealing with non-stationary data, the problem of assessing the sensitivity of trend 
inferences is particularly delicate. Beyond the serious conceptual problems explained in 
Sec. ??, one is confronted with several practical issues. Most of the sensitivity tests in 
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Part I were based on examining a shorter portion of same time series or on calculating 
the maxima on data blocks of different lengths. In the present non-stationary context, 
both operations would result in an alteration of the statistical properties of the sample 
(exactly because of the non-stationarity) and this makes comparisons somewhat ambigu- 
ous. An example is provided in Fig. ??, where we compute the best estimate GEV fits 
using sequences of yearly maxima having different lengths - but starting at the same 
instant (year 14500) - extracted from the time series with the slowest trend intensity 
IsTe = 2/ (1000 years). Notice that the best fit obtained by taking 100 yearly maxima 
is stationary. The corresponding extrapolations in time are, of course, completely wrong. 
By using 500 and 1000 maxima, the best estimates obtained (not shown) fall inside the 
confidence band of the 2000 years-based estimate for most values of time. 

The above example illustrates the trend dilemma: on the one hand, in order to be de- 
tected, a statistical trend has to be sufficiently fast with respect to the length of the record 
of observations; on the other hand, if the trend is too fast then the "adiabatic hypothesis" 
discussed in Sec. ?? is no longer valid: one is left with no "reference statistics" against 
which the inferred models can be compared. 

Moreover, when considering large time spans a further practical complication arises: 
due to the nonlinear dependence of the statistical properties with respect to the external 
parameter Te, a functional relation between the GEV parameters and time might require 
many parameters to achieve goodness-of-fit. Therefore, one faces the problem of large 
uncertainties in the parameter estimates or even lack of convergence. This has indeed 
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been observed for the present time series: if we consider a long record, such that the 
change in Te is large, the model family Gj, g with parameters as in (??) becomes inade- 
quate to catch the time dependence of the statistics of extremes. A first indication of this 
was reported at the end of Sec. ??. As a further example, we have examined a data block 
of length 5000 starting at year 14500 in the time series with AT^ = 2/(1000 years). In- 
spection of graphical diagnostics (probability and quantile plots) reveals that no model in 
the family Gp^q produces an acceptable inference. It should be emphasized that goodness- 
of-fit is achieved for the same time series using blocks of length 1000, that is, performing 
inferences that are more localized in time. So in this case the problem is not the failure of 
the "adiabatic" or smoothness hypothesis, but the nonlinear dependence of the attractor 
on the parameter Te, which manifests itself on sufficiently large time intervals. 



7. Smooth dependence on the forcing 

The set-up of the present analysis (see Sec. ??) has been chosen to allow comparison of 
the non- stationary GEV inferences with the results of Part I, obtained from statistically 
stationary time series. To perform the comparison, for each i = 1, . . . , 21 the best esti- 
mate parameters (/i*(t), a'(t), ^*) inferred from data block 5' are first expressed as func- 
tions of Te inside the interval (??). This is achieved by inverting the trend formula (??), 
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(writing time as a function of Tg): 



Te-T^+1 



Te e [9,51]. 



(21) 



and inserting this into the expression of {jl\t) ,a\t) , ^'^) . This yields functions that are 
denoted as (/i* (Tg) , a* (T^) , . These are evaluated at the central point Tg of the interval 
of definition and plotted in Fig. ??. Confidence intervals are given by the same worst 
case estimate (??) used for Fig. ??. A rather smooth dependence on Te is observed, 
especially for the GEV parameters /i and o. The location parameter /i turns out to be not 
very sensitive to changes in the trend intensity, being much more sensibly dependent on 
variations in Te- Moreover its confidence intervals are always very small (relatively to 
the size of /x). 

Denote by Go,o(iy*) the time-independent GEV models inferred in Part I from sta- 
tionary sequences of 1000 yearly maxima, computed with Te fixed at T|;. Denote as 
/i\ 5"*, and ^* the inferred values of the GEV parameters of Gqx){'w^)- Since the graphs of 
the parameters /i*, a'\ and versus very closely match those in Fig. ??, comparison 
with the stationary data is presented under the form of relative differences (Fig. ??). To 
be precise, on the left column the absolute values of the ratios 



V^{T'e)-~^' 



(22) 
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are plotted against (similarly for the GEV parameters a and Remarkable agreement 
is obtained for the parameter fi: the relative differences less than 10% and drop below 5% 
for large Te and for all considered trend intensities. Excellent agreement is also obtained 
for a (particularly for large Te) and for ^ except for the fastest trend intensity ATe = 
2/(100 years). In the latter case, indeed, the sample uncertainty is as large as (or even 
larger than) the estimates self. 

We emphasize that inferring time-independent models G'o,o('2*) from the non-stationary 
data z* would induce very large errors (particularly in the scale and shape parameters), 
compare Fig. ??. A much better (even surprising) agreement between the stationary and 
non- stationary estimates is obtained with the procedure described in the previous section: 
first fitting the time-dependent model q»(z*) and then evaluating its parameters at the 
central point T^. There is agreement even in the estimates of the parameter ^, which is 
usually the most difficult one to infer. Indeed, the inferred values are negative, in the case 
of stationary time series, since the attractor is bounded and since the energy observable 
E{t) is a continuous function of the phase space variables, the total energy is bounded on 
any orbit lying on (or converging to) the attractor. Therefore, the total energy extremes 
are necessarily WeibuU distributed (^ is negative) Part I. Although this property is not 
bound to hold for non-stationary forcing, it is still verified, see Tab. ??. 

Two distinct power law regimes are identified for the GEV parameters (/i*, a\ as 
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functions of T^, having the form 

/2'(T|) =«^(T|)> and d\Th) = a^i^V^ , (23) 

ee Fig. ?? and Fig. ??. The values of the fitted exponents 7^ and 7^ in each scaling regime 
are reported in Tab. ?? and Tab. ??, respectively. A similar power law dependence of the 
GEV parameters on Te was already observed in Part I for the stationary data sets w^: 
indeed, the exponents obtained there are very similar to those in Tab. ?? and Tab. ??, 
particularly for large Te- The lack of a power-law scaling regime for the parameter a for 
small Te explains both the more pronounced differences between the stationary and non- 
stationary estimates (Fig. ??) and the necessity of including a quadratic term in fi and/or 
a linear term for a in the statistical model to get acceptable inferences. This highlights 
the strongly nonlinear behaviour of the baroclinic model, whose response to changes of 
Te has different features depending on the considered range of variation. 

Two factors explain the qualitative analogies and the quantitative agreements between 
the time-dependent models discussed here and the stationary results of Part I. First of all, 
the trend intensity imposed on Te in all cases is sufficiently slow with respect to the time 
of relaxation of the baroclinic model to the statistics of extreme values of the total energy. 
For clarity, we emphasize that the latter time scale is that used in Sec. ?? to define the 
"adiabatic" hypothesis: it is the time necessary to obtain a good sampling of the SRB 
measure on the attractor, provided that one may consider the system as "frozen" (with 
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constant Tg) for sufficiently long time spans. We do not know whether this time scale 
bears any physical relation with other time scales, such as those of baroclinic instability 
or low-frequency variability (both have been described in ? for the present model). The 
second factor is that the system's statistical behaviour responds rather smoothly to the 
imposed time-dependent variation of the parameter Te- This smooth dependence on Te 
of the statistical properties of the baroclinic model was analyzed in detail in ?? by consid- 
ering not only global physical quantities such as total energy and average wind profiles, 
but also finer dynamical indicators, such as the Lyapunov exponents and dimension. Both 
properties of smoothness and "adiabaticity" are of crucial importance in order to justify 
the usage of non-stationary GEV models that are (locally) smooth functions of time, such 
as our polinomial family Gp,g. 



8. Summary and Conclusions 

In this paper we have proposed a general, although not universal, framework for the anal- 
ysis of trends in extremes of climatic time series. When all the shortcomings which are 
present in datasets and observations have to be considered, a rigorous definition of ex- 
tremes and a neat, clean, and legible approach to the evaluation of trends is necessary in 
order to get useful and reliable information (?). The time-dependent approach allows to 
express the inferred GEV distributional parameters as functions of time. As expected, 
it is found that trend in the statistics of extreme values is detectable in a reliable way. 
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provided that the record of observations is sufficiently long, depending on the time scale 
of the trend itself. Trend inference and assessment is much more problematic than in the 
statistically stationary inference. First of all, one is faced with a serious conceptual prob- 
lem: there is no "operational" definition of probability, since, to say it in loose words, the 
time series is not a sampling of a unique probability distribution, as it is in the stationary 
case. Even if one assumes that the time series is a realization of a sequence of random 
variables (with different distributions), the statistical properties of the sample are altered 
by any operation such as resampling or taking shorter subsamples, which makes sensitiv- 
ity studies somewhat ambiguous. One must assume that the distributions of the random 
variables vary slowly and smoothly with time, so that the time series contains sufficient 
sampling information on the "local" (in time) statistical behaviour. 

In the present context, we have adopted GEV models whose parameters are poly- 
nomial in time: the location parameter n is at most quadratic with respect to time and 
the scale parameter a is at most linear in time. Since the relation between the macro- 
scopic forcing Te and time is invertible, the time dependence of the inferred GEV mod- 
els can be expressed as a relation between the GEV parameters and Te, showing rather 
interesting properties. The location and scale parameters feature power-law dependence 
with respect to Te, while the shape parameter has in all cases a negative value. As ex- 
pected, both results are in agreement with what obtained in the companion paper (Part I) 
for stationary data. Since the parameter Te increases monotonically in the simulations 
with the baroclinic model, the system certainly does not possess any invariant measure. 
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However, the results suggest that, as Te increases, the system explores statistical states 
which vary smoothly with Te and whose properties are locally quite similar to those ob- 
tained in the stationary setting. This is even captured for the relatively fast trend intensity 
ATg = 2/(100 years). The proposed explanation is that: 

1. the system's statistical properties depend rather smoothly on Te (also compare (??)); 

2. the adopted time-scales of variation of Te {i.e., the trend intensity ATg;) are suffi- 
ciently slow compared to the relaxation time to the statistics of extreme values. 

The second condition, that was explained in more detail in Sec. ??, amounts to the heuris- 
tic statement that for sufficiently short time spans the system's statistical properties can 
be considered /rozen to those holding for a corresponding value of Te- The possibility 
of using GEV models that are locally smooth (polynomial functions of time) depends 
essentially on these two conditions. For example, for a system having several bifurca- 
tions as the control parameter is changed the time-dependent GEV modeling would be 
much more complicated. This problem is currently under investigation. However, even if 
the above two conditions do hold, the inference of time-dependent GEV models is valid 
locally in time, that is, if the sequences of maxima used for the inference span not too 
large time periods. For large time spans, indeed, the non-linear response of the baroclinic 
model to variations in Te becomes dominant and polynomial GEV models are no longer 
suitable. On the other hand, if the sequences of maxima used for the inference are too 
short (depending on the trend intensity), wrong trend estimates may be obtained. 
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We conclude by observing that the present and the companion paper (Part I) are de- 
voted not merely to the statistical inference of extremes and their trends but also to explore 
the possibility of using extreme statistics in diagnosing the dynamical state of a geophys- 
ical fluid. Our analysis of the problem reveals, in fact, that diagnostics which is based 
on "universal" (GEV theorem), robust (smoothness properties), simple (power-law scal- 
ing), controllable (low-dimensional parametric) statistical models can help very much in 
setting up well targeted models of the general circulation, see (??). 

There are several ways in which we plan to extend the present study. First of all, we 
have considered a rather global indicator, the total energy of the system. Other choices 
might be to analyze the wave kinetic energy, the available energy or also the maximum 
vorticity on the domain of the model, which might behave differently as Te is changed. 
Moreover, there are delicate issues connected with reducing the scale from a global indi- 
cator to a local one, such as the value of the wind on a grid point. This brings into play 
all complications due to the multifractality and the spatial dependence of the process. A 
further development of the present work is the usage of extreme statistics as a dynamical 
indicator, in the sense of process oriented metrics ?. All these issues are currently under 
investigation. 
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